[ <--- prev -- ] [ HOME ] [ -- next ---> ]
Unlike some other Monte Carlo particle transport codes, FLUKA gets its input
mainly from a simple file. It offers a rich choice of options for scoring most
quantities of possible interest and for applying different variance reduction
techniques, without requiring the user to write a single line of code.
However, although normally there is no need for any "user code", there are
special cases where this is unavoidable, either because of the complexity of
the problem, or because the desired information is too unusual or too
problem-specific to be offered as a standard option.
And on the other hand, even when this is not strictly necessary, experienced
programmers may like to create customised input/output interfaces.
A number of user routines (available on LINUX and UNIX platforms in directory
usermvax) allow to define non-standard input and output, and in some cases
even to modify to a limited extent the normal particle transport. Most of them
are already present in the FLUKA library as dummy or template routines, and
require a special command in the standard input file to be activated. Users
can modify any one of these routines, and even insert into them further calls
to their own private ones, or to external packages (at their own risk!). This
increased flexibility must be balanced against the advantage of using as far as
possible the FLUKA standard facilities, which are known to be reliable and well
tested.
To implement their own code, users must perform the following steps:
A typical way to do this is:
..............
LOGICAL LFIRST
SAVE LFIRST
DATA LFIRST /.TRUE./
* return message from first call
IF (LFIRST) THEN
WRITE(LUNOUT,*) 'Version xxx of Routine yyy called'
LFIRST = .FALSE.
ENDIF
..............
IMPORTANT: The user should not modify the value of any argument in a routine calling list, except when marked as "returned" in the description of the routine here below. Similarly, no variable contained in COMMON blocks should be overwritten unless explicitly indicated.
fff yyy.f (produces a new file yyy.o)
lfluka -o myfluka -m fluka yyy.o
This will produce a new executable (indicated here as myfluka). To run the new executable, launch the usual rfluka script with the option -e myfluka.
It is recommended that at least the following lines be present at the beginning of each routine:
INCLUDE '(DBLPRC)'
INCLUDE '(DIMPAR)'
INCLUDE '(IOUNIT)'
Each INCLUDE contains a COMMON block, plus related constants.
Additional INCLUDEs may be useful, in particular BEAMCM, CASLIM, EMFSTK, SOURCM,
EVTFLG, FHEAVY, GENSTK, LTCLCM, FLKMAT, RESNUC, SCOHLP, SOUEVT, FLKSTK, SUMCOU,
TRACKR, USRBIN, USRBDX, USRTRC, USRYLD.
Files flukaadd,add and emfadd.add contain a full documentation about the
meaning of the variables of these INCLUDE files.
Here is a summary of their content:
DBLPRC: included in ALL routines of FLUKA, contains the declaration
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
and sets many mathematical and physical constants. Users are strongly
encouraged to adhere to "FLUKA style" by using systematically double
precision (except for very good reasons such as calling external single
precision scoring packages), and to use constants defined in this file
for maximum accuracy.
DIMPAR: dimensions of the most important arrays
IOUNIT: logical input and output unit numbers
BEAMCM: properties of primary particles as defined by options BEAM
and BEAMPOS
CASLIM: number of primary particles followed (needed for normalisation)
EMFSTK: particle stack for electrons and photons
SOURCM: user variables and information for a user-written source
EVTFLG: event flags (still undocumented)
FHEAVY: stack of heavy secondaries created in nuclear evaporation
GENSTK: properties of each secondary created in a hadronic event
LTCLCM: LaTtice CeLl CoMmon (needed when writing symmetry transformations for
Lattice Geometry)
FLKMAT: material properties
RESNUC: properties of the current residual nucleus
SCOHLP: SCOring HeLP (information on current detector and estimator type). It
contains a flag ISCRNG, indicating the quantity being scored by the
current binning or the estimator corresponding to the current
detector, and one JSCRNG corresponding to the binning/detector number.
Binnings and detectors are sequentially numbered according to their
order of appearance in standard input. Note that detectors
corresponding to different estimators can have the same JSCRNG number
(for instance Binning N. 3 and Tracklength detector N. 3). They can be
distinguished based on the value of ISCRNG. However, note that the
same value of ISCRNG may have different meanings in functions FLUSCW
and COMSCW.
SOUEVT: SOUrce EVenT (useful when source particles are obtained from an
external event generator)
FLKSTK: main FLUKA particle stack
SUMCOU: total numbers and total weights relative to many physical and Monte
Carlo events (needed for normalisation, energy balance etc.)
TRACKR: TRACKs Recording (properties of the currently transported particle and
its path)
USRBIN, USRBDX, USRSNC, USRTRC, USRYLD: parameters of the requested estimator
detectors and binnings
Argument list (all variables are input only):
WVLNGT : photon wavelength (in cm)
OMGPHO : angular frequency (omega = 2pi nu) of the photon (in s-1)
MMAT : material index
Function ABSCFF returns a user-defined absorption coefficient for optical photons. It is activated by setting WHAT(2) < -99 in command OPT-PROP, with SDUM = blank. See option OPT-PROP and Chap. (12) for more information.
Argument list:
IJ : particle type (1 = proton, 8 = neutron, etc.: see code in (5))
Input only, cannot be modified.
XA,YA,ZA : current particle position
MREG : current geometry region
RULL : amount to be deposited (unweighted)
LLO : particle generation. Input only, cannot be modified.
ICALL : internal code calling flag (not for general use)
This routine is activated by option USERWEIG, with WHAT(6) > 0.0. Energy and star densities obtained via SCORE and USRBIN, and energy and stars obtained via EVENTBIN and production of residual nuclei obtained via RESNUCLEi are multiplied by the value returned by this function. The user can implement any desired logic to differentiate the returned value according to any information contained in the argument list (particle type, position, region, amount deposited, particle generation), or information available in the SCOHLP COMMON block (binning number, type of scored quantity). The scored quantity is given by the flag ISCRNG (in SCOHLP):
ISCRNG = 1 --> Energy density binning
ISCRNG = 2 --> Star density binning
ISCRNG = 3 --> Residual nuclei scoring
The binning/detector number is given by JSCRNG (in SCOHLP) and is printed in output:
Res. nuclei n. 3 äny-name" , "high" energy products, region n. 4
R-Phi-Z binning n. 5 öther-name" , generalised particle n. 1
Note that a detector of residual nuclei can have the same JSCRNG number as a binning (use the value of ISCRNG to discriminate). Further information can be obtained including COMMON TRACKR (for instance particle's total energy, direction cosines, age). TRACKR contains also special user variables (both integer and in double precision) which can be used to save information about particles which have undergone some particular event. If data concerning the current material are needed, it can be accessed as MEDIUM(MREG) if file (FLKMAT) is included. Indeed, a common simple application of COMSCW is to score dose according to the local density (especially useful to get the correct average dose in bins straddling a boundary between two different media):
..................
INCLUDE '(FLKMAT)'
INCLUDE '(SCOHLP)'
..................
* ========== In order to compute doses ========= *
* Medium(n) is the material number of region n
* Rho(m) is the density of material m (in g/cm3)
* Iscrng = 1 means we are depositing energy (not stars)
IF ( ISCRNG .EQ. 1 ) THEN
* to get dose in Gy (elcmks is the electron charge in C)
COMSCW = ELCMKS * 1.D12 / RHO (MEDIUM(MREG))
ELSE
* oneone is defined as 1.D0 in include DBLPRC
COMSCW = ONEONE
ENDIF
..................
Note that the variables in the argument list, with the exception of IJ, LLO and ICALL, are local copies of those used for particle transport, and therefore can be modified to have an effect on scoring, without affecting transport.
Note: setting the variable LSCZER = .TRUE. before RETURN (LSCZER is in COMMON SCOHLP), will cause zero scoring whatever the value returned by COMSCW. This is more efficient than returning a zero value.
Argument list (all variables are input only):
WVLNGT : photon wavelength (in cm)
OMGPHO : angular frequency (omega = 2pi nu) of the photon (in s-1)
MMAT : material index
Function DFFCFF returns a user-defined diffusion coefficient for optical photons. It is activated by setting WHAT(3) < -99 in command OPT-PROP, with SDUM = blank. See option OPT-PROP and Chap. (12) for more information.
Argument list (all variables are input only):
IJ : particle type (input only)
NTRUCK : number of step points (input only)
XTRUCK,YTRUCK,ZTRUCK : particle step points, can be modified by user
MREG : region number (input only)
LLO : particle generation (input only)
ICALL : internal code calling flag (not for general use)
Subroutine ENDSCP allows to shift by a user-defined distance the energy which is being deposited along a step or several step binning portions, by providing new segment endpoints. A typical application is to simulate an instrument drift.
Argument list (all variables are input only):
IJ : particle type
PLA : particle momentum (if > 0), or kinetic energy (if < 0)
(input only)
TXX,TYY,TZZ : particle direction cosines, can be modified by user
NTRUCK : number of step points
XTRUCK,YTRUCK,ZTRUCK : particle step points, can be modified by user
NREG : new region number (input only)
IOLREG : old region number (input only)
LLO : particle generation (input only)
ICALL : internal code calling flag (not for general use)
Subroutine FLDSCP allows to shift by a user-defined distance the track whose length is being scored as fluence along a step or several step binning portions, by providing new segment endpoints. A typical application is to simulate an instrument drift.
Argument list:
IJ : particle type (input only, cannot be modified)
PLA : particle momentum (if > 0), or -PLA = kinetic energy
(if < 0)
TXX,TYY,TZZ: particle current direction cosines
WEE : particle weight
XX,YY,ZZ: particle position
NREG : current region (after boundary crossing)
IOLREG : previous region (before boundary crossing). Useful only
whith boundary crossing estimators; for other estimators it
has no meaming.
LLO : particle generation (input only, cannot be modified)
ICALL : internal code calling flag (not for general use)
Function FLUSCW is activated by option USERWEIG, with WHAT(3) > 0.0. Yields obtained via USRYIELD, fluences calculated with USRBDX, USRTRACK, USRCOLL, USRBIN, and currents calculated with USRBDX are multiplied by the value returned by this function. The user can implement any desired logic to differentiate the returned value according to any information contained in the argument list (particle type, energy, direction, weight, position, region, boundary, particle generation), or information available in the SCOHLP COMMON block (binning or detector number, estimator type). The scored quantity is given by the flag ISCRNG (in SCOHLP):
ISCRNG = 1 --> Boundary crossing estimator
ISCRNG = 2 --> Track length binning
ISCRNG = 3 --> Track length estimator
ISCRNG = 4 --> Collision density estimator
ISCRNG = 5 --> Yield estimator
The binning/detector number is given by JSCRNG (in SCOHLP) and is printed in output:
Bdrx n. 2 "bdxname" , generalised particle n. 8,
from region n. 22 to region n. 78
Track n. 6 "trkname" , generalised particle n. 14, region n. 9.
Note that a track-length detector can have the same JSCRNG number as
a boundary crossing one or a binning etc. (use the value of ISCRNG to
discriminate).
Further information can be obtained including COMMON TRACKR (for
instance particle age). TRACKR contains also special user variables
(both integer and in double precision) which can be used to save
information about particles which have undergone some particular event.
Function FLUSCW has many applications. A common one in shielding
calculations is to multiply selected scored fluences by
particle/energy-dependent fluence-to-dose equivalent conversion
coefficients, or by some instrument response, radiation damage curve,
etc. A version of FLUSCW converting particle fluences to various
forms of dose equivalent and effective dose is available in file
deq99.f, in the $FLUPRO/flutil directory. Instructions for its use
are contained in the comment lines at the beginning of the file.
Another application is conditional scoring (score only if
within a certain distance from a point, etc.): for instance it is
possible to implement a sort of 2-dimensional fluence binning on
a plane boundary.
FLUSCW can be used also when scoring "fluxes" (i.e. fluences or
currents or yields) of heavy ions, to discriminate between different
types of ions.
All ions in FLUKA carry the same id-number IJ = -2: therefore ion
fluxes obtained by the various estimators will refer to all ions
unless an additional discrimination is introduced by the user at
scoring time. This can be done in FLUSCW as follows:
CALL USRDCI(IJ,IONA,IONZ,IONM)
The three integer values returned are the following ion properties:
IONA = mass number of the ion
IONZ = atomic number
IONM = flag for isomeric state
Based on their values, the user can decide or not to accept the
scoring.
Other interesting applications are based on the fact that FLUSCW is called at every boundary crossing, provided that at least one USRBDX detector has been requested. Although the function has been designed mainly to weight scored quantities, it can be "cheated" to do all sorts of side things, even not directly connected with scoring. Note that the variables in the argument list, with the exception of IJ, LLO and ICALL, are local copies of those used for particle transport, and therefore can be modified to have an effect on scoring, without affecting transport.
Note: setting the variable LSCZER (in SCOHLP) = .TRUE. before
RETURN, will cause zero scoring whatever the value returned
by FLUSCW. This is more efficient than returning a zero value.
Argument list (all variables are input only):
IJ : particle code, except that it is set to 3 for both e+ and e-
QU2 : squared momentum transfer (GeV/c)^2
ZMEDIU : atomic number of target nucleus
AMEDIU : atomic mass of target nucleus
Function FORMFU can be used to override the standard value of the nuclear charge form factor. It must return the squared value of the nuclear charge form factor for particle IJ. The default version computes the form factor in Born approximation for a medium of given composition, using the simple expression given by Tsai [Tsa74], and accounts also for the contribution of incoherent scattering. The function is called by the multiple and single scattering routines if option MULSOPT has been issued with WHAT(3) <t 0.0 for electrons and positrons, and WHAT(2) < 0.0 for hadrons and muons. See Note 2) to option MULSOPT.
Argument list (all variables are input only):
TXX, TYY, TZZ : particle direction cosines
UXSRFC, UYSRFC, UZSRFC: direction of the normal to the surface
MREG : region the particle is coming from
NEWREG : region the particle is going to
MMAT : material the particle is coming from
MMATNW : material the particle is going to
Function FRGHNS can be used to return a non-zero value for the roughness of a boundary between two materials, relevant for optical photon transport (default roughness is zero for all boundaries). Meaningful only if options OPT-PROP or OPT-PROD have been requested.
These three functions are used to define 3-dimensional fluence distributions to be calculated by special user-defined binnings (see USRBIN with WHAT(1) = 8.0 in the first card).
Argument list (all variables are input only):
IJ : particle type
PCONTR : particle momentum
XA,YA,ZA : particle position
MREG : current region
LATCLL : current lattice cell
ICALL : internal code calling flag (not for general use)
Argument list: same as for MUSRBR above
Argument list: same as for MUSRBR above (except LATCLL)
The 3 functions are called at tracklength events. What is scored is
the particle track-length multiplied by the particle's weight,
possibly modified by a user-written FLUSCW (see above).
Subroutine LATTIC is activated by one or more LATTICE cards in the
geometry input (see (8)). It is expected to transform coordinates and
direction cosines from any lattice cell (defined by card LATTICE) to
the reference system in which the basic structure has been described.
The user is expected to provide a transformation of coordinates and
vector direction cosines from each lattice cell to the corresponding
basic structure (in ENTRY LATTIC) and of direction cosines from the
basic structure to each corresponding lattice cell (in ENTRY LATNOR).
Entries:
Argument list:
XB(1), XB(2), XB(3) : actual physical position coordinates in IRLTGG
lattice cell
WB(1), WB(2), WB(3) : actual physical direction cosines in IRLTGG
lattice cell
DIST : current step length
SB(1), SB(2), SB(3) : transformed coordinates in prototype cell
UB(1), UB(2), UB(3) : transformed cosines in prototype cell
IR : region number in prototype cell
IRLTGG : lattice cell number
IRLT : array containing region indices corresponding to
lattice cells
IFLAG : reserved variable
LATTIC returns the tracking point coordinates (SB) and direction cosines(UB) in the reference prototype geometrical structure, corresponding to real position/direction XB, WB in the actual cell IRLTGG (defined as input region IR by a LATTICE card). When the lattice option is activated, the tracking proceeds in two different systems: the "real" one, and that of the basic symmetry unit. Particle positions and directions are swapped from their real values to their symmetric ones in the basic cell, to perform the physical transport in the regions and materials that form the prototype geometrical structure and back again to the real world. The correspondence between "real" and "basic" position/direction depends on the symmetry transformation and on the lattice cell number.
Argument list:
UN(1), UN(2), UN(3) : direction cosines of the vector normal to the
surface, in the prototype cell (entry values)
and in the lattice cell (returned values)
IRLTNO : present lattice cell number
Entry LATNOR transforms the direction cosines stored in the vector
UN(3) from the system of the basic prototype unit to that of the real
world in lattice cell number IRLTNO. Therefore, this cosine
transformation must be the inverse of that performed on the cosines by
the LATTIC entry: but while LATTIC maps vector UB to a different
vector WB, LATNOR maps the UN vector to itself.
Note that if the transformation implies a rotation, it is necessary to
save first the incoming UN cosines to local variables, to avoid
overwriting the vector before all transformation statements are executed.
Notes:
Different symmetry transformations can of course be implemented
in the same LATTIC routine (each being activated by a different cell
number or range of cell numbers).
The advantage of the lattice geometry is to avoid describing in detail
the geometry of repetitive multi-modular structures. It must be
realised, however, that a penalty is generally paid in computer efficiency.
Also, a region contained in the prototype cell and all those "mapped"
to it inside lattice cells are treated by the program as if they
were connected with "non-overlapping ORs" into a single region.
Therefore, any region-based scoring (options SCORE, USRTRACK, etc.)
can only provide quantities averaged over the whole structure. More
detailed information must be obtained by region-independent options
such as USRBIN or by user-written routines (MGDRAW).
The USRBIN and EVENTBIN options, with WHAT(1) = 8, can also be used
to request a special binning type which activates the MUSRBR,
LUSRBL, FUSRBV user routines to recover lattice information (see
routine musrbr.f below)
A transformation between a lattice cell and a prototype region can
alternatively be defined without resorting to the LATTIC user routine.
In this case, the transformation is defined via a ROT-DEFIni card and
the correspondence is established by giving the transformation index
in the SDUM of the LATTICE card (see Chap. (8)).
Argument list:
X,Y,Z : current position (input only)
BTX,BTY,BTZ : direction cosines of the magnetic field vector
(returned).
B : magnetic field intensity in Tesla (returned)
NREG : current region (input only)
IDISC : if returned = 1, the particle will be discarded
MAGFLD is activated by option MGNFIELD with WHAT(4-6) = 0.0 and is used to return intensity and direction of a magnetic field based on the current position and region. It is called only if the current region has been flagged as having a non-zero magnetic field by option ASSIGNMAt, with WHAT(5) = 1.0. The magnetic field spatial distribution is often read and interpolated from an external field map. Note that in any case the direction cosines MUST be properly normalised in double precision (e.g. BTX = SQRT(ONEONE - BTY**2 - BTZ**2)), even if B = 0.0. Please read carefully the notes on option MGNFIELD.
Argument list:
IFLAG : type of nuclear interaction which has produced secondaries:
1 : inelastic
2 : elastic
3 : low-energy neutron
NUMSEC : number of secondary particles produced in the interaction
MDSTCK is called after a nuclear interaction in which at least one secondary particle has been produced, before any biasing is applied to decide which secondary will be loaded in the main stack for further transport. The properties of the secondaries are stored in the secondary stack (COMMON GENSTK). With MDSTCK, the user can analyse those secondaries, write them to a file, or even modify the content of GENSTK (for instance applying his own biasing). In the latter case, however, it is his responsibility to make sure that energy is conserved, the various physical quantities are still consistent, etc.
Subroutine MGDRAW, activated by option USERDUMP with WHAT(1) >= 100.,
usually writes a "collision tape", i.e. a file where all or selected
transport events are recorded. The default version (unmodified by the
user) offers several possibilities, selected by WHAT(3) in USERDUMP.
Details are given in (11).
Additional flexibility is offered by a user entry USDRAW, selected
by WHAT(4) in USERDUMP, interfaced with the most important physical
events happening during particle transport.
The user can modify of course also any other entry of this
subroutine (BXDRAW called at boundary crossings, EEDRAW called at event
end, MGDRAW for trajectory drawing, ENDRAW for recording of
energy depositions and SODRAW for recording of source events:
for instance the format of the output file can be changed, and
different combinations of events can be written to file.
No information is written by default at EEDRAW and BXDRAW calls, but the
entries are called for any value of WHAT(3) in USERDUMP (EEDRAW also
for WHAT(4) >= 1).
But the most interesting aspect of the routine is that the five
entries (all of which, if desired, can be activated at the same time
by setting USERDUMP with WHAT(3) = 0.0 and WHAT(4) >= 1.0) constitute
a complete interface to the whole FLUKA transport. Therefore, MGDRAW
can be used not only to write a collision tape, but to do any kind
of complex analysis (for instance studying correlations between events).
Entries:
Argument list (all variables are input only):
ICODE : FLUKA physical compartment originating the call
= 1: call from Kaskad (hadrons and muons)
= 2: call from Emfsco (electrons, positrons and photons)
= 3: call from Kasneu (low-energy neutrons)
= 4: call from Kashea (heavy ions)
= 5: call from Kasoph (optical photons)
MREG : current region
MGDRAW writes by default, for each trajectory, the following variables (contained in COMMON TRACKR):
NTRACK : number of track segments
MTRACK : number of energy deposition events along the track
JTRACK : type of particle
ETRACK : total energy of the particle
WTRACK : weight of the particle
Ntrack values of XTRACK, YTRACK, ZTRACK : end of each track segment
Mtrack values of DTRACK : energy deposited at each deposition event
CTRACK : total length of the curved path
Other variables are available in TRACKR (but not written by MGDRAW unless the latter is modified by the user: particle momentum, direction cosines, cosines of the polarisation vector, age, generation, etc. (see a full list in the comment in the INCLUDE file).
Argument list (all variables are input only):
ICODE : FLUKA physical compartment originating the call
= 19: call from Kaskad (hadrons and muons)
= 29: call from Emfsco (electrons, positrons and photons)
= 39: call from Kasneu (low-energy neutrons)
= 49: call from Kashea (heavy ions)
= 59: call from Kasoph (optical photons)
MREG : number of region before boundary crossing
NEWREG : number of region after boundary crossing
XSCO, YSCO, ZSCO : coordinates of crossing point
BXDRAW is called at each boundary crossing (if requested by the user
with USERDUMP, WHAT(3) < 7.0). There is no default output: any
output must be supplied by the user.
If name-based input is being used, the names corresponding to MREG
and NEWREG can be obtained via a call to routine GEOR2N:
CALL GEOR2N (NUMREG, NAMREG, IERR)
where NUMREG (input variable) is a region number, and NAMREG (returned variable) is the corresponding region name (to be declared as CHARACTER*8). IERR is a returned error code: if = 0, the conversion is successful. Example:
.......................................
CHARACTER*8 MRGNAM, NRGNAM
.......................................
ENTRY BXDRAW ( ICODE, MREG, NEWREG, XSCO, YSCO, ZSCO )
CALL GEOR2N ( MREG, MRGNAM, IERR1 )
CALL GEOR2N ( NEWREG, NRGNAM, IERR2 )
IF(IERR1 .NE. 0 .OR. IERR2 .NE. 0) STOP "Error in name conversion"
.......................................
IF(MRGNAM .EQ. "MyUpsREG" .AND. NRGNAM .EQ. "MyDwnREG") THEN
.......................................
Argument list:
ICODE = -1: event not completed
= 0: normal event termination
= 4: stack overflow
EEDRAW is called at the end of each event, or primary history, (if requested by the user with USERDUMP, WHAT(3) =< 0.0). There is no default output: any output must be supplied by the user.
Argument list (all variables are input only):
ICODE : type of event originating energy deposition
1x: call from Kaskad (hadrons and muons)
10: elastic interaction recoil
11: inelastic interaction recoil
12: stopping particle
13: pseudo-neutron deposition
14: particle escaping (energy deposited in blackhole)
15: time kill
2x: call from Emfsco (electrons, positrons and photons)
20: local energy deposition (i.e. photoelectric)
21 or 22: particle below threshold
23: particle escaping (energy deposited in blackhole)
24: time kill
3x: call from Kasneu (low-energy neutrons)
30: target recoil
31: neutron below threshold
32: neutron escaping (energy deposited in blackhole)
33: time kill
4x: call from Kashea (heavy ions)
40: ion escaping (energy deposited in blackhole)
41: time kill
42: delta ray stack overflow
5x: call from Kasoph (optical photons)
50: optical photon absorption
51: optical photon escaping (energy deposited in
blackhole)
52: time kill
MREG : current region
RULL : energy amount deposited
XSCO, YSCO, ZSCO : point where energy is deposited
ENDRAW writes by default, for each energy deposition point:
0 : flag identifying ENDRAW output from that of other entries
ICODE : see argument list
JTRACK, ETRACK, WTRACK : see MGDRAW above
XSCO, YSCO, ZSCO, RULL : see argument list
No arguments
SODRAW writes by default, for each source or beam particle:
-NCASE (in COMMON CASLIM, with a minus sign to identify
Sodraw output): number of primaries followed so far
NPFLKA (in COMMON FLKSTK): stack pointer
NSTMAX (in COMMON FLKSTK): highest value of the stack pointer
encountered so far
TKESUM (in COMMON SOURCM): total kinetic energy of the
primaries of a user written source (see source.f
here below), if applicable. Otherwise = 0.0
WEIPRI (in COMMON SUMCOU): total weight of the primaries
handled so far
Npflka times: ILOFLK : type of source particle
TKEFLK+AM : total particle energy (kinetic+mass)
WTFLK : source particle weight
XFLK, YFLK, ZFLK : source particle position
TXFLK, TYFLK, TZFLK : source particle direction cosines
Argument list (all variables are input only):
ICODE : 10x: call from Kaskad (hadron and muon interactions)
100: elastic interaction secondaries
101: inelastic interaction secondaries
102: particle decay secondaries
103: delta ray generation secondaries
104: pair production secondaries
105: bremsstrahlung secondaries
110: radioactive decay products
20x: call from Emfsco (electron, positron and photon
interactions)
208: bremsstrahlung secondaries
210: Moller secondaries
212: Bhabha secondaries
214: in-flight annihilation secondaries
215: annihilation at rest secondaries
217: pair production secondaries
219: Compton scattering secondaries
221: photoelectric secondaries
225: Rayleigh scattering secondaries
30x: call from Kasneu (low-energy neutron interactions)
300: interaction secondaries
40x: call from Kashea (heavy ion interactions)
400: delta ray generation secondaries
MREG : current region
XSCO, YSCO, ZSCO : interaction point
USDRAW is called after each particle interaction (if requested by
the user with USERDUMP, WHAT(4) >= 1.0). There is no default output:
any output must be supplied by the user.
Information about the secondary particles produced is available
in COMMON GENSTK, except that concerning delta rays produced by heavy
ions (in which case the properties of the single electron produced
are available in COMMON EMFSTK, with index NP).
Another exception is that about heavy evaporation fragments (deuterons,
3-H, 3-He, alphas, with JTRACK ID equal respectively to -3, -4, -5, -6)
and fission/fragmentation products generated in an inelastic interaction
(with JTRACK = -7 to -12) which are all stored in COMMON FHEAVY, with
index NPHEAV. To get the kinetic energy of particles with JTRACK < -6,
one must subtract from their total energy (ETRACK} in COMMON TRACKR)
their fully stripped nuclear mass (AMNHEA in COMMON FHEAVY).
Information about the interacting particle and its trajectory can be
found in COMMON TRACKR (see description under the MGDRAW entry
above). In TRACKR there are also some spare variables at the user's
disposal: LLOUSE (integer), ISPUSR (integer array) and SPAUSR
(double precision array). Like many other TRACKR variables, each of
them has a correspondent in the particle stacks, i.e. the COMMONs
from which the particles are unloaded at the beginning of their
transport: FLKSTK, EMFSTK and OPPHST (respectively, the stack of
hadrons/muons, electrons/photons, and optical photons). The
correspondence with TRACKR is shown below under STUPRF/STUPRE.
When a particle is generated, its properties (weight, momentum,
energy, coordinates etc., as well as the values of the user flags)
are loaded into one of the stacks. The user can write a STUPRF or
STUPRE subroutine (see description below) to change anyone of such
flags just before it is saved in stack.
When a particle starts to be transported, its stack variables are
copied to the corresponding TRACKR ones. Unlike the other TRACKR
variables, which in general become modified during transport due to
energy loss, scattering etc., the user flags keep their original
value copied from stack until they are changed by the user himself
(generally in USDRAW).
One common application is the following: after an interaction which
has produced sencondaries, let USDRAW copy some properties of the
interacting particle into the TRACKR user variables. When STUPRF is
called next to load the secondaries into stack, by default it copies
the TRACKR user variables to the stack ones. In this way,
information about the parent can be still carried by its daughters
(and possibly by further descendants). This technique is sometimes
referred to as "latching".
Argument list (all variables are input only):
OMGPHO : angular frequency (omega = 2 pi nu) of the photon (in s-1)
WVLNGT : photon wavelength (in cm)
MREG : old region number
NEWREG : new region number
SIGANW : absorption coefficient in the new region (cm-1)
SIGDNW : diffusion coefficient in the new region (cm-1)
RFNDPR : refractive index in the new region
VGRPNW : group velocity in the new region (cm s-1)
LPHKLL : if .TRUE., the photon will be absorbed on the boundary
Subroutine OPHBDX sets the optical properties of a boundary surface. The call is activated by command OPT-PROP, with SDUM = SPEC-BDX. See option OPT-PROP and Chap. (12) for more information.
Argument list (all variables are input only):
WVLNGT : photon wavelength (in cm)
OMGPHO : angular frequency (omega = 2pi nu) of the photon (in s-1)
Function QUEFFC returns a user-defined quantum efficiency for an optical photon of the given wavelength or frequency. It is activated by OPT-PROP, with SDUM = SENSITIV, by setting the 0th photon sensitivity parameter to a value < -99. See option OPT-PROP and Chap. (12) for more information.
Argument list (all variables are input only):
WVLNGT : photon wavelength (in cm)
OMGPHO : angular frequency (omega = 2pi nu) of the photon (in s-1)
MMAT : material index
Function RFLCTV returns a user-defined value equal to 1-R, where R is the reflectivity of the current material for an optical photon of the given wavelength or frequency. It is activated by OPT-PROP, with SDUM = METAL, and WHAT(3) < -99. See option OPT-PROP and Chap. (12) for more information.
Argument list (all variables are input only):
WVLNGT : photon wavelength (in cm)
OMGPHO : angular frequency (omega = 2pi nu) of the photon (in s-1)
MMAT : material index
Function RFRNDX returns a user-defined refraction index of the current material for an optical photon of the given wavelength or frequency. It is activated by OPT-PROP, with SDUM = blank, and WHAT(1) < -99. See option OPT-PROP and Chap. (12) for more information.
No arguments
Subroutine SOEVSV is always called after a beam particle is loaded
into FLKSTK, but a call to SOEVSV can be inserted by the user anywhere
in a user routine.
SOEVSV copies the whole FLKSTK to another COMMON, SOUEVT, which can
be included in other user routines. In other words, this routine is
used to "take a snapshot" of the particle bank at a particular time
for further use (interfacing to independent generators, etc.)
Argument list:
NOMORE : if set = 1, no more calls will occur (the run will be
terminated after exhausting the primary particles loaded
into FLKSTK stack in the present call). The history number
limit set with option START will be overridden.
Subroutine SOURCE is probably the most frequently used user routine. It is activated by option SOURCE and is used to sample primary particle properties from distributions (in space, energy, time, direction, polarisation or mixture of particles) too complicated to be described with the BEAM, BEAMPOS and POLARIZAti cards alone. For each phase-space variable, a value must be loaded into COMMON FLKSTK (particle bank) before returning control. These values can be read from a file, generated by some sampling algorithm, or just assigned.
Reading from a file is needed, for instance, when the particle data are taken from a collision file, written by FLUKA or by another program. The user must open the file with a unit number > 20 (unit numbers lower than 20 are reserved), in one of the following ways:
Then, a READ statement in SOURCE can be used to get the data to load in stack, for instance:
READ(21,*) IPART, X, Y, Z, COSX, COSY, COSZ, ENERGY, WEIGHT
ILOFLK (NPFLKA) = IPART
XFLK (NPFLKA) = X
YFLK (NPFLKA) = Y
ZFLK (NPFLKA) = Z
TXFLK (NPFLKA) = COSX
...etc...
(NPFLKA is the current stack index).
Direct assignment can be done explicitly, for instance:
PMOFLK (NPFLKA) = 305.2D0
or implicitly, leaving unmodified values input with BEAM or BEAMPOS:
PMOFLK (NPFLKA) = PBEAM
(PBEAM is the momentum value input as WHAT(1) in option BEAM). A set of direct assignments, one for each of several different stack entries, can be useful, for example, to define a series of RAYs through the geometry (see (14)):
DO 10 I = 1, 20
NPFLKA = NPFLKA + 1
ILOFLK (NPFLKA) = 0 (0 is the RAY particle id number)
XFLK (NPFLKA) = 500.D0 + DBLE(I) * 40.D0
YFLK (NPFLKA) = 200.D0
...etc...
10 CONTINUE
To sample from a uniform distribution, the user must use the function FLRNDM(DUMMY), which returns a double precision pseudo-random number uniformly distributed between 0 (included) and 1 (not included). Actually, DUMMY can be any variable name. A simple example of sampling from a uniform distribution is that of a linear source along the Z axis, between Z = 10 and Z = 80:
Z1 = 10.D0
Z2 = 80.D0
ZFLK (NPFLKA) = 10.D0 + (Z2 - Z1) * FLRNDM(XXX)
One way to sample a value XX from a generic distribution f(x) is the following. First integrate the distribution function, analytically or numerically, and normalise to 1 to obtain the normalised cumulative distribution:
/x /xmax
F(x) = | f(x)dx / | f(x)dx
/xmin /xmin
Then, sample a uniform pseudo-random number t using FLRNDM and get the desired result by finding the inverse value:
-1
XX = F (t)
(analytically or most often by interpolation).
A FLUKA subroutine is available to sample directly from a Gaussian
distribution:
CALL FLNRRN (RGAUSS)
or, if two independent Gaussian distributed numbers are needed:
CALL FLNRR2 (RGAUS1, RGAUS2)
(faster than calling FLNRRN twice).
The technique for sampling from a generic distribution described above can be extended to modify the probability of sampling in different parts of the interval (importance sampling). We replace f(x) by a weighted function g(x) = f(x) * h(x), where h(x) is any appropriate function of x we like to choose. We normalise g(x) in the same way as f(x) before:
/x /xmax /x
G(x) = | g(x)dx / | g(x)dx = | f(x)*h(x)dx / B
/xmin /xmin /xmin
and we need also the integral of f(x) over the whole interval:
/xmax
A = | f(x)dx
/xmin
All the sampling is done using the biased cumulative normalised function G instead of the original unbiased F: we sample a uniform pseudo-random number t as before, and we get the sampled value XX by inverting G(x):
-1
XX = G (t)
The particle is assigned a weight B/(A * h(XX))
A special case of importance sampling is when the biasing function
chosen is the inverse of the unbiased distribution function:
h(x) = 1/f(x)
g(x) = f(x) * h(x) = 1
G(x) = (x - xmin) / (xmax - xmin)
In this case we sample a uniform pseudo-random number t using FLRNDM as shown above. The sampled value XX is simply given by:
XX = xmin + (xmax - xmin)*t
and the particle is assigned a weight:
/xmax
B/(A h(X)) = f(XX)*(xmax - xmin) / | f(x)dx
/xmin
But since FLUKA normalizes all results per unit primary weight, any
constant factor is eliminated in the normalization. Therefore it is
sufficient to assign each particle a weight f(x).
Because XX is sampled with the same probability over all possible
values of x, independently of the value f(XX) of the function,
this technique is used to ensure that sampling is done uniformly
over the whole interval, even though f(x) might have very small
values somewhere. For instance it may be important to avoid
undersampling in the high-energy tail of a spectrum, steeply falling
with energy but more penetrating, such as that of cosmic rays or
synchrotron radiation.
Option SOURCE allows the user to input up to 12 numerical values
(WHASOU(1),(2)...(12)) and one 8-character string (SDUSOU) which can
be accessed by the subroutine by including the following line:
INCLUDE '(SOURCM)'
These values can be used as parameters or switches for a multi-source
routine capable to handle several cases, or to identify an external
file to be read, etc., without having to compile and link again the
routine.
In the SOURCE routine there are a number of mandatory statements,
(clearly marked as such in accompanying comments) which must not be
removed or modified. The following IF block initialises the total
kinetic energy of the primary particles and sets two flags: the
first to skip the IF block in all next calls, and the second to
remind the program, when writing the final output, that a user source
has been used:
* +